/*
  File:    AngiogneesisInterface.cpp
  Purpose: This provides abstracted interfaces to the angiogenesis module.

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <cmath>
#include <ctime>

using namespace std;

#include "AngiogenesisInterface.h"

#include "LogfileEntry.h"
#include "ParseText.h"
#include "Matrix.h"
#include "ImprovedMath.h"

#include "UCI_TumorGrowth/Oxygen.h"
#include "UCI_TumorGrowth/Pressure.h"

Matrix Pressure_Previous;
VesselStructure MyVessel;
CellVesselStructure TheVessel;

// from driver


extern char* StringBuffer;
extern char* LogFile;
extern char* ParameterFile;

bool ForceBloodflowWithoutLooping;

// put any globals necessary from the angio here

// from angiogenesis
extern double AngioDeltaT;
extern double AngioGrowthInterval;
extern double AngioTime;


// from tumor growth

extern double dt;
extern double dx;
extern int Xnodes;
extern int Ynodes;
extern double t;
extern Matrix TAF;

// for triggering
extern double AngiogenesisTriggeringThreshold;
bool AngiogenesisIsStarted = false;
bool AngiogenesisDisabled = false;
bool IfPutArtery = false;
bool AngiogenesisResume = false;
extern bool b_venousIncluded;
extern bool FlowBegin;
extern bool AdaptionBegin;

extern bool TopBloodVessel;
extern bool BottomBloodVessel;

double pressure_extravasation( double tumor_pressure , double vessel_pressure )
{
 extern double collapseConstant;

 extern double Custom_Cutoff_Max;
// Custom_Cutoff_Max = 1.0;//0.7;//RetrieveDoubleParameter( ParameterFile , "Custom_Cutoff_Max" ); 
 
// collapseConstant = 3375.0;//RetrieveDoubleParameter( ParameterFile , "collapseConstant" ); 
 
 double dp = ( collapseConstant*tumor_pressure - vessel_pressure )/(Custom_Cutoff_Max*collapseConstant);
 
 if( dp < 0.0 )
 { dp = 0.0; }
 if( dp > 1.0 )
 { dp = 1.0; }
 
 
// double output = ( 3.0*IntPow( input , 2) - 2.0*IntPow( input , 3 ) ); // smooth cubic

 double output = dp; // linear 
 
 // At this point, output is the cutoff factor. We want extravasation: 
 
 return 1.0-output; 
}

bool InitializeAngiogenesis( void )
{
 if( AngiogenesisDisabled ) 
 { return false; }

 if( AngiogenesisIsStarted == true )
 {
  StartStringBuffer( "Warning! Angiogenesis already initialized. Ignoring request!\n" );
  cout << StringBuffer;
  MakeLogEntry( LogFile, StringBuffer );
  return false;
 }
 
 extern double BloodTissueTransferRate;
 extern double OxygenUptakeRate;
 extern double Oxygen_pi;

 Oxygen_pi = RetrieveDoubleParameter( ParameterFile , "Oxygen_pi" );
 BloodTissueTransferRate = RetrieveDoubleParameter( ParameterFile , "BloodTissueTransferRate" );
 AngiogenesisTriggeringThreshold = RetrieveDoubleParameter( ParameterFile , "AngiogenesisTriggeringThreshold" );
 OxygenUptakeRate =  RetrieveDoubleParameter( ParameterFile , "OxygenUptakeRate" );
// cout<<OxygenUptakeRate;
 //getchar();
 ForceBloodflowWithoutLooping = (bool) RetrieveDoubleParameter( ParameterFile, "ForceBloodflowWithoutLooping" );
 
 // put any external function calls to initialize anigogenesis here
 
 // MIN'S OnRun_AngioInitialization();
 
TheVessel.InitializeCellTime(AngiogenesisResume);

TheVessel.InitializeCellList(AngiogenesisResume);
//cout<<"pass"<<endl;
//getchar();
CellVesselStructure* Vessel;
Vessel = TheVessel.GetAddress();

MyVessel.VesselInitialization( AngiogenesisResume, Vessel );




  AngiogenesisIsStarted = true;
  return true;
}

bool CheckForAngiogenesisTrigger( void )
{
 if( AngiogenesisDisabled ) 
 { return false; }

 if( AngiogenesisIsStarted == true )
 { return false; }
 
 // Check to see if we should resume:
 // If the simulation is resuming and angio resume data exists, trigger the 
 // standard initialization routine and overwrite with resume data. 
 // Otherwise, continue with the normal trigger check. 
 
 extern bool ResumeFromPreviousData;
 if( ResumeFromPreviousData )
 {
  extern char SaveDirectory [1025];
  char resume_file_to_check [1024];
  sprintf( resume_file_to_check , "%s/angio_resume.dat" , SaveDirectory );

  FILE* angio_temp_file = fopen( resume_file_to_check , "r" );
  if( angio_temp_file != NULL )
  { 
   AngiogenesisResume = true;
   fclose( angio_temp_file );
   return true;  
  }
 }
 
 AngiogenesisTriggeringThreshold = RetrieveDoubleParameter( ParameterFile , "AngiogenesisTriggeringThreshold" );
 // double Tolerance = square( dx*AngiogenesisTriggeringThreshold );
 double Tolerance = AngiogenesisTriggeringThreshold;
 
 // check the gradients along the parent vessel
 
 /* Fix to look at value instead of gradient */
 /*******
 for( int i=1 ; i < Xnodes-1 ; i++ )
 {
  /*
  double dGradX = 0; // 0.5*( *TAF(i+1,Ynodes-1) - *TAF(i-1,Ynodes-1) );
  double dGradY = ( *TAF(i,Ynodes-1) - *TAF(i,Ynodes-2) );
  double dTemp = square( dGradX ) + square( dGradY );
     */
	 /****
  double dTempTop = *TAF(i,Ynodes-1);
  double dTempBottom = *TAF(i,0);
  
  if( !BottomBloodVessel )
  { dTempBottom = -1.0; }
  if( !TopBloodVessel )
  { dTempTop = -1; }
  
  double dTemp = MAX( dTempTop, dTempBottom);
  
  if( dTemp >= Tolerance )
  {
   StartStringBuffer( "Angiogenesis triggered at time t = " );
   sprintf( StringBuffer + strlen(StringBuffer) , "%3.8f\n" , t );
   cout << StringBuffer;
   MakeLogEntry( LogFile, StringBuffer );
   return true; 
  }
 }
 ******/
 
{
	StartStringBuffer( "Angiogenesis triggered at time t = " );
    sprintf( StringBuffer + strlen(StringBuffer) , "%3.8f\n" , t );
    cout << StringBuffer;
    MakeLogEntry( LogFile, StringBuffer );
    return true; 

} 

 
 return false;
}



bool UpdateAngiogenesis( void )
{
 extern Matrix BloodVesselIndicator;
 extern Matrix BloodVesselIndicatorX;
 extern Matrix BloodVesselIndicatorY;
 extern Matrix HaematocritLevel;
 extern Matrix BloodVesselRadius;
 extern Matrix BloodNodePressure;
 extern Matrix Pressure;
 extern Matrix Flow;
 extern Matrix SproutTips;
 extern Matrix SproutAge;
 extern Matrix TAF;
 extern Matrix MDE;
 extern Matrix ECM;
	extern Matrix Shuffle;
	//extern double t;
 extern bool ResumeFromPreviousData;
 extern bool AdaptionBegin;
 extern bool ResetPressure;

 //cout << "UpdateAngiogenesis start, flowbegin: " << FlowBegin << endl;
 //getchar();

 if( AngiogenesisDisabled ) 
 { return false; }

 //Min's Notes 1: put the Artery Vessel.
 //Also we can change the code here to change the
 //Position or shape of the Atery.
 //Question is if we are not putting Artery at
 //Boundary, we should change the BC of the O2 solver.
 
 
 if( !IfPutArtery && !ResumeFromPreviousData)
 {
	
	//Just give position, other factors will be defined later on.
	//If want to change position, change function ArteryPosition
	MyVessel.ArteryPosition();
	
	IfPutArtery = 1;
	
 }

 if( AngiogenesisIsStarted == false )
 {
  if( !CheckForAngiogenesisTrigger() )
  { return false; }
  else
  {
   InitializeAngiogenesis();
  }
 }


 if( BloodVesselIndicator.TellRows() != Xnodes || 
     BloodVesselIndicator.TellCols() != Ynodes )
 { BloodVesselIndicator = zeros(Xnodes,Ynodes); }
 if( BloodVesselIndicatorX.TellRows() != Xnodes || 
     BloodVesselIndicatorX.TellCols() != Ynodes )
 { BloodVesselIndicatorX = zeros(Xnodes,Ynodes); }
 if( BloodVesselIndicatorY.TellRows() != Xnodes || 
     BloodVesselIndicatorY.TellCols() != Ynodes )
 { BloodVesselIndicatorY = zeros(Xnodes,Ynodes); }
 if( BloodVesselRadius.TellRows() != Xnodes || 
     BloodVesselRadius.TellCols() != Ynodes )
 { BloodVesselRadius = zeros(Xnodes,Ynodes); }
 
 if( Flow.TellRows() != Xnodes || 
     Flow.TellCols() != Ynodes )
 { Flow = zeros(Xnodes,Ynodes); }
 if( SproutTips.TellRows() != Xnodes || 
     SproutTips.TellCols() != Ynodes )
 { SproutTips = zeros(Xnodes,Ynodes); }
 if( SproutAge.TellRows() != Xnodes || 
     SproutAge.TellCols() != Ynodes )
 { SproutAge = zeros(Xnodes,Ynodes); }
 



 Shuffle = zeros(Xnodes,Ynodes);
  extern double Custom_Cutoff_Max_New;
  extern double collapseConstant_New;

 Custom_Cutoff_Max_New = RetrieveDoubleParameter( ParameterFile , "Custom_Cutoff_Max_New" ); 
 if( Custom_Cutoff_Max_New < 1e-6 )
	{ Custom_Cutoff_Max_New = 1.0; }


 collapseConstant_New =  RetrieveDoubleParameter( ParameterFile , "collapseConstant_New" ); 
 //ResetPressure = 0;
 
 for(int i=0;i<Xnodes;i++)
 {
	for(int j=0;j<Ynodes;j++)
	{
	
	//if(i>80&&i<120&&j>80&&j<120&&((t>0.05&&t<0.15)||(t>4&&t<4.1)||(t>4.5&&t<4.6)||(t>6&&t<6.1)||(t>3.5&&t<3.6)))
	//{
	//*Shuffle(i,j)=1;
	//	cout<<i<<"   "<<j<<endl;
	//}
	
	
		
	//if(( collapseConstant_New* *Pressure(i,j) - *BloodNodePressure(i,j) )/(Custom_Cutoff_Max_New*collapseConstant_New)>1.0)
	if(*Pressure(i,j) > collapseConstant_New)
	{
	*Shuffle(i,j)= 1;
	//ResetPressure = 1;
	}
	
	
	}
 }
 
 //cout<<"here"<<endl;
//getchar();
CellVesselStructure* Vessel;
Vessel = TheVessel.GetAddress();



  MyVessel.UpdateVesselSystemSteadyStateShuffle(Vessel);
  //UpdateRadiusByPressure();
  MyVessel.SaveVesselList();// MIN'S OnRun_WriteAngioVessels();

// TheVessel.GrowSprout(TAF);
 //TheVessel.String_Interact(TAF);
 //TheVessel.GenerateMatrix();
 //TheVessel.SaveVesselList();// MIN'S OnRun_WriteAngioVessels();


 
 
 return true;
}

double TransferOxygenFromBloodVessels( double o2, double p, int i, int j )
{ 
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix Phi;
 extern double _EPS;
 extern double dx;
 extern double dy;
 extern double BloodTissueTransferRate;
 extern Matrix BloodVesselIndicator;
 extern bool TumorBegin;
 extern bool FlowBegin;
 extern Matrix HaematocritLevel;
 extern Matrix Flow;
 double HaematocritScale = 0.5; 
 double HaematocritToOxygenScale = 10;
 extern Matrix VesselCells;
 double FlowRef=1.9096e-11; //m3/s (Flow in PV)
 double FlowAverage = 2.8e-16;
 extern Matrix BloodNodePressure;

 extern Matrix BloodVesselPressureX;
 extern Matrix BloodVesselPressureY;
 if( HaematocritLevel.TellRows() != Xnodes ||
     HaematocritLevel.TellCols() != Ynodes )
 { HaematocritLevel = zeros(Xnodes,Ynodes); }
/* 
double f = *Flow(i,j)/FlowAverage;
 if(f>1)
 {
 f=1;
 }

 
 double h = *HaematocritLevel(i,j)* *Flow(i,j);
*/ 
 
 double h = *HaematocritLevel(i,j)/0.45;

 
 extern Matrix Pressure;
/*
 double extravasation = pressure_extravasation( *Pressure(i,j) , *BloodVesselPressureX(i,j) );

 if( *BloodVesselIndicator(i,j) > 0.25 && ForceBloodflowWithoutLooping == true )
 { h = 1.0; extravasation = 1.0; }

 if( h < 0.05 ) 
 { return 0.0; }

*/
double extravasation = pressure_extravasation( *Pressure(i,j) , *BloodNodePressure(i,j) );
h = h * extravasation;
if(*BloodVesselIndicator(i,j)>0.25)
{
if(!FlowBegin)
{
return BloodTissueTransferRate;
}
else
{
 //return f*BloodTissueTransferRate;
 return h*BloodTissueTransferRate;
} 
 
 
}
else
{
return 0;
}
}

void UpdateRadiusByPressure()
{
    double RadiusLimitInf=RetrieveDoubleParameter( ParameterFile , "RadiusLimitInf" ); 
	double RadiusLimitSup=RetrieveDoubleParameter( ParameterFile , "RadiusLimitSup" );
 extern Matrix BloodVesselRadiusXdir;
 extern Matrix BloodVesselRadiusYdir;
 extern Matrix BloodVesselRadius;
 extern Matrix BloodVesselIndicator;
 extern Matrix Pressure;
 extern Matrix BloodVesselPressureX;
  extern Matrix BloodVesselPressureY;
 for(int i=0; i<Xnodes; i++)
 {
	for(int j=0; j<Ynodes; j++)
	{
	if(*BloodVesselIndicator(i,j)>0)
	{
	*BloodVesselRadiusXdir(i,j)=RadiusLimitInf+(*BloodVesselRadiusXdir(i,j)-RadiusLimitInf)*pressure_extravasation( *Pressure(i,j), *BloodVesselPressureX(i,j) );
	*BloodVesselRadiusYdir(i,j)=RadiusLimitInf+(*BloodVesselRadiusYdir(i,j)-RadiusLimitInf)*pressure_extravasation( *Pressure(i,j), *BloodVesselPressureY(i,j) );
	*BloodVesselRadius(i,j)=RadiusLimitInf+(*BloodVesselRadius(i,j)-RadiusLimitInf)*pressure_extravasation( *Pressure(i,j), *BloodVesselPressureX(i,j) );
	}
	}
 }


}

